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Abstract 

Monte Carlo simulations of 2D vortex lattice melting in a thin superconduct- 
ing film (or alternatively an array of Josephson junctions) are performed in 
the London limit. Finite size scaling analyses are used to make a detailed 
test of the dislocation mediated melting theory of KTNHY. We find that the 
melting transition is weakly first order, with a jump in shear modulus very 
close to that predicted by the KTNHY theory. No hexatic liquid phase is 
found. 

PACS number: 64.60-i, 74.60-w, 74.76-w 
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Interest in the melting transition of two dimensional (2D) vortex lattices has revived 
recently, primarily due to the belief that the strongly fluctuating, layered, high-T^ super- 
conductors may behave 2D like in sufficiently large magnetic fields This 2D melting 
transition has generally been believed to be described by the Kosterlitz-Thouless-Nelson- 
Halper in- Young (KTNHY) theory [^|-|^ of dislocation mediated melting. However, the very 
existence of a 2D vortex lattice at any finite temperature has been recently questioned by 
Moore p, who first argued that phase fiuctuations destroy such a lattice, and then found 
support for this picture from Monte Carlo (MC) simulations P]. High order, high tempera- 
ture perturbative expansions [0], similarly show no evidence for freezing into a vortex lattice 
in 2D. Other MC simulations P,^, however, find clear evidence for a melting transition at 
finite temperature. Hu and MacDonald [§] find this transition to be first order, in opposition 
to the KTNHY prediction. 

The above cited simulations, have all been performed in the "lowest Landau level" ap- 
proximation, in which the complex order parameter ip{r) is expanded in terms of eigenstates 
of the Gaussian part of the Landau-Ginsburg free-energy functional. Alternatively, one may 
use instead the London approximation, in which the amplitude of ipi^) is assumed to be 
constant, and only the phase is allowed to vary. In this limit, the problem can be efficiently 
simulated by utilizing the well known mapping |lO| onto the 2D Coulomb gas. Logarithmi- 



cally interacting point charges model vortices in the phase of ip{r). For a uniform magnetic 
field B, one has a fixed density -B/$o of positive integer charges, on a uniform neutralizing 
background ($o is the fiux quantum). The London approximation should be valid whenever 
the bare vortex core radius (^v, the average spacing between vortices. This corresponds 
to temperatures well below the mean field transition temperature, where vortex lattice melt- 
ing is expected to occur. Earlier simulations of this 2D Coulomb gas, in the context of the 



2D one component plasma problem show clear evidence for a finite temperature melt- 
ing transition, and suggest that the transition is weakly first order. In the following, we 
report on new simulations in this London approximation, in which we carry out the first 
finite size scaling analysis of the melting transition, making a detailed comparison with the 



KTNHY theory. We show that the shear modulus jumps discontinuously to zero at the 
melting transition with a value very close to the universal KTNHY prediction, however we 
find no evidence for a hexatic liquid phase. We perform the first conclusive test of the order 
of the melting transition (within the London approximation) using the histogram method. 
We find that the transition is weakly first order, consistent with earlier suggestions. 
The model we simulate is given by the Hamiltonian, 

^ = ^E(^^-/)(%-/)nr.-r,). (1) 

In order to reduce the size of phase space, we have discretized the continuum by constraining 
the charges to the sites i of a periodic triangular grid with spacing oq. The sum is over all 
pairs of sites in an L x L parallelogram-shaped cluster of triangular grid, and Ui = 0, 
or +1, is the point charge (vorticity) at site i. The neutralizing background charge is 
/ = {\/3al/2)(B /(^o) = (ao/a^,)^, and V{r) is the lattice Coulomb potential in 2D, which 
solves, 

AV(r) = -27r5r,o (2) 

subject to periodic boundary conditions. Here is the discrete Laplacian. To keep the 
energy finite, it is necessary to preserve total charge neutrality, which leads to the constraint 
Nc = I]i = where iV = is the total number of sites in the grid, and is the total 
number of charges in the system. Thus / is the density of charges. Further details may 
be found in Ref. ||12|. The connection between Eq. ([^) and the superconducting system is 



obtained by measuring the temperature of the Coulomb gas model, in units of ^qcZ/Stt^A^, 
where A is the magnetic penetration length, and d the thickness, of the superconducting 
layer 

The Hamiltonian (|I]) has been studied extensively [l^ll^ in the dense charge limit (/ = 



1/2, 1/3). Here we are interested in the dilute limit / <^ 1, where we expect our discretized 
model to well approximate the continuum (see Refs. |]TT| for similar simulations, directly in 



the continuum). We always choose / commensurate with the system length L, so that the 



ground state will be a perfect triangular charge (vortex) lattice. We study various charge 
densities / = with m = 3 to 12, and fixed Nc ~ 100. Detailed finite size scaling 

analyses are carried out for the specific case of / = 1/49, and Nc = 16, 25, ... , 169. 

Our MC updating scheme is as follows. In each MC step one charge is selected at random 
and moved to a different site within a radius ~ This excitation is then accepted or 

rejected according to the standard Metropolis algorithm. Nc such attempts we refer to as 
one MC sweep. At low temperature, we also make global moves, by attempting to shift 
entire rows of charges by one space. Such moves are meant to model long wavelength shear 
excitations, and help to accelerate equilibration near the vortex lattice melting transition. 
Data is collected by heating the system up from the ground state. At each temperature we 
discard 30, 000 MC sweeps to equilibrate the system. Then, starting from this equilibrated 
configuration, we perform several (typically 4 — 6) independent runs of 100, 000 sweeps each 
to sample physical quantities. Errors are estimated from the standard deviation of these 
independent runs. To verify consistency of our results, we also perform cooling from a 
random configuration at high temperature; no substantial hysteresis is found. 

The physical quantities we measure are: (i) the inverse dielectric function. 



where rik = Si exp(— ik ■ r^). The vanishing of e ^ upon heating signals an "insulator- 
conductor" transition in the Coulomb gas. As e^^ can be mapped onto the helicity modulus 



In the simulation, the k — limit is approximated by averaging e ^ over the three smallest 
allowed wave vectors; (u) the six-fold orientational order correlation. 



where the sum is over sites with non- vanishing charges rii = +1, and 6i is the angle of the 
bond from Ui to its nearest neighbor, relative to some fixed reference direction; and {in) the 
structure function 




(3) 



of the superconductor [0], its vanishing signals the loss of superconducting phase coherence. 




(4) 
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In Fig. 1, we plot e^^(T) and feiT), versus T, for / = 1/49 and Nc = 169. The behavior 
of ipe{T) indicates two separate transitions at Tc{f) and Tm- For a simple visualization of 
the resulting three phases, we show in Fig. 2 intensity plots of 5'(k), for k's in the first 
Brillouin zone (BZ). We show results for three different temperatures, with the data for 
each value of T restricted to one third of the BZ. For T = 0.003 (Fig. 2a), just below Tc{f), 
we see a regular array of (5-function Bragg peaks, indicating long range translational order. 
Thus for T < Tc{f), we have a vortex lattice which is pinned to the discretizing grid. For 
T = 0.0065 (Fig. 2b), just below Tm, we see a regular array of peaks, but the peaks are 
now of finite width. These peaks are consistent with power law singularities, characteristic 
of the algebraic translational order expected for a 2D lattice in the continuum. Thus for 
Tcif) <T < we have a "floating" vortex lattice, which is depinned from the grid, and 
we have reached the continuum limit. For T = 0.0075 (Fig. 2c), slightly above T^, we see a 
rotationally invariant structure (y^g ~ 0), typical for a liquid with short range correlations. 
Thus for T > T^, the floating lattice has melted into a liquid. 

Returning to Fig. 1, we see that e'^ vanishes at the depinning transition Tc{f). Thus 
the floating lattice has lost superconducting phase coherence. This is nothing more than a 
reflection of the flux flow resistance to be expected from an unpinned vortex lattice, which 
is free to drift transversly to an applied d.c. current. Our results explicitly show that the 
absence of phase coherence in this k ^ sense, does not not imply the absence of a well 
defined vortex lattice. 

In the inset to Fig. 1, we show the dependence of Tc{f) and Tm on the charge density /. 
We see that only for sufficiently dilute systems, /< 1/25, is there a floating lattice phase; for 
/ > 1/25 there is only a single transition from a pinned lattice to a liquid. As / decreases, 
Tcif) — * as ~ /, while Tm quickly approaches a flnite constant Tm = 0.0070 ± 0.0005, in 
good agreement with the melting temperature found in earlier continuum simulations |TT[]. In 
terms of the superconductor, this means a vortex lattice melting at Tm = 0.0070 ^q'^/Stt^A^, 
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well within the bounds estimated by Fisher [Q] from the KTNHY theory. 

The transition at Tc(/) is an artifact of our discretization of the continuum, and hence 
has no direct physical meaning for a uniform continuous superconductor. However Tc(/) 
does represent a physical depinning transition for the related problem of vortex states in 
periodic superconducting networks, such as Josephson junction arrays. In this case, our 
result that Tc{f —* 0) vanishes, is consistent with early commensurability arguments by 
Teitel and Jayaprakash The result that a floating lattice exists above Tc(/) is, however, 
a new observation in the Josephson array context; the melting of this lattice at a finite 
may dominate the physics of such arrays at small /. 

To investigate the nature of the melting transition T^, we have carried out detailed 
finite size scaling analyses for the case / = 1/49, in which is well separated from 
Tc. Our approach is guided by the KTNHY theory 0. For a 2D lattice in the contin- 
uum, translational correlations decay algebraically with a temperature dependent expo- 
nent, (e**^'('''^'"j'') ~ |ri — rj\^''T^^'^\ where G is a reciprocal lattice vector of the real space 
charge lattice. For the 2D superconducting case, where the vortex compressibility is infinite, 
Vg{T) = k bT IGl"^ /Att^, where fi is the vortex shear modulus. If Gi is the shortest reciprocal 
lattice vector, then the KTNHY theory predicts that at T^, 77gi takes a discontinuous jump 
to infinity from the universal value of riG-^(T^) = 1/3. 

To test this prediction for translational order, we measure the height of peaks in the 
structure function. From Eq. (^) we find that these should scale as 

^(G) ~ L^-^^C^) for T<Tm. (6) 

Above Tm, translational order has exponential decay with a correlation length ^. One then 
obtains 

S{G)^e for T>T^. (7) 

In the Fig. 3a we plot S{Gi)/ L"^ , as a function of L on a log-log scale, for several different 
temperatures. Data for each temperature fall on a straight line, confirming the expected 
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power-law behavior. These straight hnes fall into three distinct groups. For T <Tc — 0.0045, 
S{Gi)/L'^ ~ 1, indicating the long range order of the pinned lattice. For T^. < T < — 
0.007, we find algebraic decay, S{Gi)/L^ ~ L~'>^^^\ For T > T^, we find S{Gi)/L'^ ~ L'^ 
with X ^ 2 as T increases, consistent with the short range order of a liquid. The lines in Fig. 
3a are a fit to Eq.(^); the resulting exponents riG-^^(T) are shown in Table 1. We see that ?7gi 
first exceeds the KTNHY universal value of 1/3 at T = 0.0065, very close to our estimated 
melting transition of Tm — 0.007, where the slopes of the lines in Fig. 3a show an apparent 
discontinuous jump. Similar results, within the "lowest Landau level" approximation, have 
very recently been obtained by Sasik and Stroud [|^. 

As a consistency check, we have also computed S'(G2), where G2 = 2Gi. Using similar 
fits as in Fig. 3a, we determine the exponent rjc^, and show the results in Table. 1. We see 
that ?7g2 — 4r7Gi as expected, since ?7g ~ |Gp. 

We now consider the orientational order. Below Tm, KTNHY predict long range 6— fold 
orientational order given by (e*®'-^*^'"-'^^*^*^-'-') ~ ae~^^^^ + 9?^. For <ti L, one obtains from 
Eq.® 

cp,r^2na(^)\ip^. (8) 



L. 

Above Tm, KTNHY predict a hexatic liquid phase, with algebraic orientational order 
^e«6(9(r)-e(o))^^ ^ y^-m(T) ^j^-j^ ?76(T) < 1/4. In such a case, one would have, 

^^^L-^^^ (9) 

At higher temperature, KTNHY predict a transition from the hexatic to an isotropic liquid, 
with short ranged orientational order. In this case, Eq.@ again holds, but with ip'^ = 



| r?| . In Fig. 3b we display (Pq{T) as a function of system size L for various temperatures. 
In the fioating solid phase below one clearly sees saturation of ^Pq{T) to a finite value as 
L increases. Solid lines represent a least square fits to Eq. (|), and the extracted values of 
ip^ are shown in Table 1. Above Tm, we try fits to both Eqs. and (Q). We always find 
that Eq. (|^) gives the superior fit, and for T > Tm we find, within our numerical precision, 
< = 0. 



Thus our results for the floating lattice phase are consistent with expectations for a 2D 
continuum lattice, and we flnd that translational correlations at the melting transition are 
consistent with the KTNHY prediction. However we do not find any evidence for a hexatic 
liquid above T^. 

The absence of the hexatic liquid suggests the possibility that the melting transition is 
not of the KTNHY type, but is perhaps weakly first order as found by Hu and MacDonald 
IP], and as suggested in Refs. lTl]|. To examine this possibility, we measured the energy 
distribution P{E) ~ e~^^^^/^ at the melting temperature Tm [|l8l, and in Fig. 4 we plot the 
resulting free energy F{E) versus E. We see a double well structure with an energy barrier 
AF between two coexisting phases. In the inset to Fig. 4 we plot the dependence of AF 
on system length L. The growth in AF as L increases is a clear signal that the transition 
is first order, although our sizes remain too small, and our data too noisy, to see clearly 
the predicted scahng AF ~ L. To determine the distributions in Fig. 4, we have computed 
P{E) at fixed T ~ T^,, and then extapolated ||T9| to determine P{E) at nearby T, finding 
the precise value of T that gives equal minima in F{E). In this way we obtain an improved 
estimate — 0.0066. 

To conclude, our results demonstrate that there is a clear finite temperature melting 
transition of the vortex lattice in two dimensions, within the London approximation. This 
transition is first order, with melting directly into an isotropic vortex liquid; no hexatic 
liquid is found. The first order transition however is very weak, so that the jump at melting 
in rjG (and hence in the vortex lattice shear modulus) remains very close to the KTNHY 
universal prediction. 

The authors are grateful to T. Chen, D. A. Huse, D. R. Nelson, and Z. Tesanovic for 
useful discussions. This work was supported by DOE grant DE-FG02-89ER14017. One of us 
(M.F.) acknowledges the Rush Rhees Fellowship of the University of Rochester for support 
in the initial stages of this project. 
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T 








0.00475 


0.188 ±0.008 


0.704 ±0.055 


0.571 ±0.007 


0.00500 


0.207 ±0.007 


0.806 ±0.032 


0.529 ±0.005 


0.00525 


0.211 ±0.007 


0.852 ±0.028 


0.504 ±0.004 


0.00550 


0.248 ±0.005 


0.998 ±0.019 


0.476 ±0.003 


0.00575 


0.255 ±0.008 


0.999 ±0.029 


0.458 ±0.003 


0.00600 


0.296 ±0.006 


1.065 ±0.028 


0.426 ±0.007 


0.00625 


0.319 ±0.010 


1.191 ±0.016 


0.403 ±0.004 


0.00650 


0.4 ±0.16 


1.4 ±0.22 


0.33 ±0.030 


0.00675 


1.4 ±0.31 


2.0 ±0.31 


0.20 ±0.041 


0.00750 


3.4 ±0.37 


3.4 ±0.44 


0.03 ±0.046 


0.01100 


2.8 ±0.23 


2.9 ±0.30 


-0.01 ±0.032 


0.01500 


2.2 ±0.12 


2.1 ±0.22 


0.00 ±0.020 



Table 1: Temperature dependence of the exponents rjG^T) and rja^iT). Also displayed 
limiting values of (p6{T) for L — cxd, ip"^. 
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FIGURES 



Fig.l Inverse dielectric function e^^(T) and orientational order correlation ^PqIT) versus T 
for / = 1/49 and Nc = 169. Inset shows the dependence of the depinning and melting 
temperatures, Tc and Tm, on charge density /. Solid and dashed lines are guides to 
the eye only. 

Fig.2 Structure function S'(k) in the first Brillouin zone, (BZ) for / = 1/49 and A^^^ = 63, 
and three different temperatures T. Data for each T is restricted to one third of the 
BZ. (a) T = 0.003, just below Tc, in the "pinned lattice" state, (b) T = 0.0065, just 
below Tm, in the "floating lattice" state, (c) T = 0.0075, just above Tm, in the liquid. 
Intensities are plotted nonlinearly to enhance features. 

Fig.3 (a) Finite size scaling of S{Gi)/L'^ (note the log-log scale). Sohd and dashed lines 
are fits to Eqs. (|^). (b) Finite size scaling of (pe{T). Solid and dashed lines are fits to 
eq. d). 

Fig.4 Free energy distribution F{E) versus E, at melting T^, for several system sizes L. 
The growth in energy barrier AF with increasing L (see inset) indicates a first order 
transition. Curves for different L are offset from each other by a constant, for the sake 
of clarity. 
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